Measuring degradation on schroeders
Fine scale acoustic perception in zebra finches
Source code and data found at https://github.com/maRce10/acoustic-fine-features-zebra-finch
1 Purpose
- Evaluate if schroeder structure is degraded during transmission
- Evaluate differencial degradation in open and closed environments
2 Report overview
3 Automatic annotation
Use the package baRulho to fin the position of schroeders in the re-recorded files
Code
master_annotations <- imp_raven(path = "./data/processed", files = "schroeder_master.wav.txt", all.data = TRUE, warbler.format = TRUE)
master_annotations$sound.id <- paste0("f0:", master_annotations$f0, "_comp:", master_annotations$label)
master_annotations$sound.id[1] <- "start_marker"
master_annotations$sound.id[nrow(master_annotations)] <- "end_marker"
names(master_annotations)
# table(master_annotations$sound.id)
markers_in_tests <-
find_markers(X = master_annotations, # annotations of sounds in master file
path = path,
) # path to supplementary files
aligned_tests <-
align_test_files(X = master_annotations, # annotations of sounds in master file
Y = markers_in_tests, # position of markers in test files
path = path, # folder with files
output = "data.frame")
check_sels(aligned_tests, path = path)
aligned_tests$start[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] <- aligned_tests$start[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] - 1.8
aligned_tests$end[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] <- aligned_tests$end[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] - 1.8
exp_raven(X = aligned_tests, file.name = "pl", path = path, sound.file.path = path)
aligned_tests <- aligned_tests[grep("marker", aligned_tests$sound.id, invert = TRUE), ]
aligned_tests_est <- selection_table(aligned_tests,extended = TRUE, confirm.extended = FALSE, path = path)
saveRDS(aligned_tests_est, file.path(path, "extended_sel_table_degradation_exp.RDS"))4 Measuring schroeder dissimilarity
4.1 Dynamic-time warping pairwise distance
- Both Schroeders have the same length
- One is duplicated and the other one is slide across the duplicated one
- The minimum DTW distance is kept as a dissimilarity measure
Code
mean_segment <- function(wave, cores = 1, plot = TRUE, pb = TRUE,
thinning = 1, col = wave_col, mean = TRUE, type = "ac", npeak = 20) {
# thin
if (thinning < 1) {
if (length(wave@left) * thinning < 10) {
stop2("thinning is too high, no enough samples left for at least 1 sound file")
}
# reduce size of envelope
wavefrm <- stats::approx(x = seq(0, duration(wave), length.out = length(wave@left)),
y = wave@left, n = round(length(wave@left) * thinning),
method = "linear")$y
} else {
wavefrm <- wave@left
}
# get empirical mode decomposition
if (type == "EMD") {
emds <- EMD::emd(wavefrm, seq_len(length(wavefrm)), boundary = "wave")
perd <- emds$imf[, 4]/max(emds$imf[, 4])
# plot(x = seq_len(length(wavefrm)), y = perd, type =
# 'l') lines(y = wavefrm / max(wavefrm), x =
# seq_len(length(wavefrm)), col = 'gray', lty = 2)
}
if (type == "ac") {
ac <- acf(x = wavefrm, lag.max = length(wavefrm), type = "covariance",
demean = FALSE, plot = FALSE)
perd <- ac$acf/max(ac$acf)
}
tpks <- seewave::fpeaks(cbind(seq_len(length(perd)), perd), plot = FALSE,
threshold = 0.5)
if (nrow(tpks) > npeak) {
tpks <- tpks[1:npeak, ]
}
segment_df <- data.frame(selec = seq_len(nrow(tpks)), pos = tpks[,
1], peak = tpks[, 2])
# get mean number of sample between peaks
mean_dist_peak <- round(mean(diff(segment_df$pos)))
segment_df$start <- segment_df$pos - mean_dist_peak/2
segment_df$end <- segment_df$pos + mean_dist_peak/2
# fix if values are out of wavefrm size
if (segment_df$start[1] > 0) {
segment_df$start[1] <- 0
}
if (segment_df$end[nrow(segment_df)] > length(wavefrm)) {
segment_df$end[nrow(segment_df)] <- length(wavefrm)
}
# extract segments into a list
segments <- lapply(seq_len(nrow(segment_df)), function(x) {
wavefrm[segment_df$start[x]:segment_df$end[x]]
})
# make all the same number of samples
segments <- lapply(segments, function(x) {
approx(x, n = max(sapply(segments, length)))$y
})
# normalize between 1, -1
segments <- lapply(segments, function(x) {
x/max(x)
})
# put all segments in a data frame
segments <- as.data.frame(segments, col.names = seq_len(length(segments)))
# compute mean segment
mean_segment <- rowMeans(segments)
if (plot) {
mean_segment_df <- data.frame(time = seq(0, 1, length.out = nrow(segments)),
mean.amp = rowMeans(segments), sd.amp = apply(segments,
1, sd))
gg <- ggplot(data = mean_segment_df, mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 25)
print(gg)
}
if (mean) {
return(mean_segment)
} else {
return(segments)
}
}Code
est_schr <- readRDS(file.path(path, "extended_sel_table_degradation_exp.RDS"))
mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
function(x) {
wave <- read_wave(est_schr, index = x)
seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
type = "ac", thinning = 0.8))
return(seg)
})
names(mean_schroeders) <- paste0(ifelse(grepl("inside", est_schr$sound.files),
"in_", "out_"), est_schr$sound.id)
saveRDS(mean_schroeders, file.path(path, "mean_schroeders_degradation_experiment.RDS"))Code
mean_schroeders <- readRDS(file.path(path, "mean_schroeders_degradation_experiment.RDS"))
mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]
mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
sd.amp = apply(mean_schroeders[[x]], 1, sd))
})
mean_schroeders_df <- do.call(rbind, mean_schroeders_list)4.2 Mean schroeders inside playback
- mean period +/- standard deviation using autocorrelation
- first 40 schroeders are shown
Code
mean_schroeders_in_df <- mean_schroeders_df[grep("in_", mean_schroeders_df$schroeder),
]
ggplot(data = mean_schroeders_in_df[mean_schroeders_in_df$schroeder %in%
unique(mean_schroeders_in_df$schroeder)[1:40], ], mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
facet_wrap("~ schroeder", ncol = 5, scales = "free_y")
Code
mean_schroeders_in <- mean_schroeders[grep("in_", names(mean_schroeders))]
nms <- names(mean_schroeders_in)
cmbs <- t(combn(nms, 2))
min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
s1 <- rowMeans(mean_schroeders_in[[cmbs[x, 1]]])
s2 <- rowMeans(mean_schroeders_in[[cmbs[x, 2]]])
# make same length if (length(s1) != length(s2))
s1 <- approx(s1, n = 100)$y
s2 <- approx(s2, n = 100)$y
# duplicate 1
s1 <- rep(s1, 2)
# run dtw over longer vector
dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]
dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
return(dtw_dist[1, 2])
}, FUN.VALUE = numeric(1))
return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})
min_dists <- do.call(rbind, min_dist_l)
min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))
names(min_dists) <- c("schr1", "schr2", "dist")
min_dists$dist <- as.numeric(min_dists$dist)
saveRDS(min_dists, file.path(path, "dtw_distance_schroeders_degradation_experiment_inside.RDS"))4.3 Mean schroeders outside playback
Code
mean_schroeders_out <- mean_schroeders_df[grep("out_", mean_schroeders_df$schroeder),
]
ggplot(data = mean_schroeders_out[mean_schroeders_out$schroeder %in%
unique(mean_schroeders_out$schroeder)[1:40], ], mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
facet_wrap("~ schroeder", ncol = 5, scales = "free_y")
Code
mean_schroeders_out <- mean_schroeders[grep("out_", names(mean_schroeders))]
nms <- names(mean_schroeders_out)
cmbs <- t(combn(nms, 2))
min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
s1 <- rowMeans(mean_schroeders_out[[cmbs[x, 1]]])
s2 <- rowMeans(mean_schroeders_out[[cmbs[x, 2]]])
# make same length if (length(s1) != length(s2))
s1 <- approx(s1, n = 100)$y
s2 <- approx(s2, n = 100)$y
# duplicate 1
s1 <- rep(s1, 2)
# run dtw over longer vector
dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]
dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
return(dtw_dist[1, 2])
}, FUN.VALUE = numeric(1))
return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})
min_dists <- do.call(rbind, min_dist_l)
min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))
names(min_dists) <- c("schr1", "schr2", "dist")
min_dists$dist <- as.numeric(min_dists$dist)
saveRDS(min_dists, file.path(path, "dtw_distance_schroeders_degradation_experiment_outside.RDS"))5 Statistical analysis
Modeling: - Multiple Regression on distance Matrices - Model:
\[\begin{align*}
Dissimilarity &\sim frequency + components + sign
\end{align*}\] - Response values scaled to make effect sizes comparable across models - Predictors were coded as pairwise binary matrices in which 0 means that calls in a dyad belong to the same level and 1 means calls belong to different levels - One model for each environment treatment as well as on the original model sounds
5.1 Original (model) Schroeders
5.2 Outside playback
Code
min_dists <- readRDS(file.path(path, "dtw_distance_schroeders_degradation_experiment_outside.RDS"))
min_dists$dist <- scale(min_dists$dist)
dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)
freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"_"), "[[", 2)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
"", sapply(strsplit(rownames(dist_tri), "_"), "[[", 3))))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"_"), "[[", 4)))
rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")
dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_outside_experiment.RDS")Code
$coef
dtw_dist pval
Int -0.017062308 0.0001
frequency 0.010099084 0.0001
components 0.001829854 0.6014
sign 0.013210074 0.0001
$r.squared
R2 pval
5.547013e-05 1.000000e-04
$F.test
F F.pval
0.7641249 0.0001000
$n
[1] 2885.3 Inside playback
Code
min_dists <- readRDS(file.path(path, "dtw_distance_schroeders_degradation_experiment_inside.RDS"))
min_dists$dist <- scale(min_dists$dist)
dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)
freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"_"), "[[", 2)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
"", sapply(strsplit(rownames(dist_tri), "_"), "[[", 3))))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"_"), "[[", 4)))
rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")
dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_inside_experiment.RDS")Code
$coef
dtw_dist pval
Int -1.0128809 1e-04
frequency 0.3789239 1e-04
components 0.2454399 1e-04
sign 0.9017830 1e-04
$r.squared
R2 pval
0.2208956 0.0001000
$F.test
F F.pval
4070.1899 0.0001
$n
[1] 2945.4 Combined results
Code
mods <- list(dtw_wv_mod = dtw_wv_mod, dtw_wv_in_mod = dtw_wv_in_mod,
dtw_wv_out_mod = dtw_wv_out_mod)
names(mods) <- c("Original Schroeders", "Inside playback", "Outside playback")
estimates <- do.call(rbind, lapply(seq_along(mods), function(x) {
Y <- data.frame(mods[[x]]$coef[-1, ])
Y$rel_coef <- Y[, 1]/max(Y[, 1])
Y$mod <- names(mods)[x]
Y$predictor <- rownames(Y)
names(Y) <- c("coef", "p", "rel_coef", "model", "predictor")
return(Y)
}))
estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
0)
estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")
estimates$model <- factor(estimates$model, levels = c("Outside playback",
"Inside playback", "Original Schroeders"))
ggplot(estimates, aes(x = predictor, y = model, fill = rel_coef)) +
geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
high = viridis(10)[7], guide = "none") + geom_text(aes(label = round(coef,
3), color = signif)) + scale_color_manual(values = c("black",
"gray")) + labs(x = "", y = "", color = "P value") + theme_classic() +
theme(axis.text.x = element_text(color = "black", size = 11, angle = 30,
vjust = 0.8, hjust = 0.8))
6 Takeaways
- Schroeder structure fairly detectable at ~30 cm
- Schroeder discrimination mostly affected in outside playback
- Harmonic content seems to be the most affect feature, both inside and outside
Session information
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] numform_0.7.0 ecodist_2.0.9 PhenotypeSpace_0.1.0
[4] ggplot2_3.4.2 baRulho_1.1.0 ohun_1.0.0
[7] Rraven_1.0.13 viridis_0.6.3 viridisLite_0.4.2
[10] warbleR_1.1.28 NatureSounds_1.0.4 tuneR_1.4.4
[13] seewave_2.2.1 formatR_1.11 knitr_1.43
[16] kableExtra_1.3.4
loaded via a namespace (and not attached):
[1] colorspace_2.1-0 rjson_0.2.21 deldir_0.2-10
[4] class_7.3-22 xaringanExtra_0.7.0 rstudioapi_0.14
[7] proxy_0.4-27 spatstat.data_2.1-0 farver_2.1.1
[10] Deriv_4.1.3 remotes_2.4.2 fansi_1.0.4
[13] xml2_1.3.5 codetools_0.2-18 splines_4.1.0
[16] polyclip_1.10-0 jsonlite_1.8.7 vdiffr_1.0.5
[19] packrat_0.9.0 cluster_2.1.2 png_0.1-7
[22] spatstat.sparse_2.0-0 compiler_4.1.0 httr_1.4.4
[25] backports_1.4.1 assertthat_0.2.1 Matrix_1.5-4.1
[28] fastmap_1.1.1 cli_3.6.1 htmltools_0.5.5
[31] tools_4.1.0 igraph_1.5.0.1 gtable_0.3.3
[34] glue_1.6.2 dplyr_1.0.10 Rcpp_1.0.11
[37] raster_3.4-13 vctrs_0.6.3 svglite_2.1.0
[40] nlme_3.1-162 xfun_0.39 stringr_1.5.0
[43] brio_1.1.3 testthat_3.1.10 rvest_1.0.3
[46] lifecycle_1.0.3 goftest_1.2-2 MASS_7.3-60
[49] scales_1.2.1 spatstat.core_2.3-0 spatstat.utils_2.2-0
[52] parallel_4.1.0 Sim.DiffProc_4.8 yaml_2.3.7
[55] pbapply_1.7-2 gridExtra_2.3 rpart_4.1-15
[58] stringi_1.7.12 e1071_1.7-13 checkmate_2.2.0
[61] permute_0.9-5 rlang_1.1.1 pkgconfig_2.0.3
[64] systemfonts_1.0.4 dtw_1.23-1 bitops_1.0-7
[67] evaluate_0.21 lattice_0.21-8 tensor_1.5
[70] sf_1.0-14 labeling_0.4.2 htmlwidgets_1.5.4
[73] tidyselect_1.2.0 magrittr_2.0.3 R6_2.5.1
[76] fftw_1.0-7 generics_0.1.3 DBI_1.1.3
[79] pillar_1.9.0 withr_2.5.0 mgcv_1.8-42
[82] units_0.8-2 abind_1.4-5 RCurl_1.98-1.12
[85] sp_1.5-1 tibble_3.2.1 crayon_1.5.2
[88] KernSmooth_2.23-21 utf8_1.2.3 spatstat.geom_2.2-2
[91] rmarkdown_2.23 grid_4.1.0 sketchy_1.0.2
[94] vegan_2.5-7 digest_0.6.33 classInt_0.4-9
[97] webshot_0.5.4 signal_0.7-7 munsell_0.5.0